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ABSTRACT 

Reverberation mapping aims to use time-delayed variations in photoionised emis- 
sion lines to map the geometry and kinematics of emission-line gas in the vicinity of an 
active galactic nucleus. By fitting to a time-variable emission line profile, one can re- 
construct a 2-dimcnsional map ^(t, v), where r is the time delay and v is the Doppler 
shift, for each emission line. In this paper we develop quasar tomography, which com- 
bines the time delay and velocity information with photoionisation physics in order to 
map the reprocessing region using information assembled from many different emis- 
sion lines. The observed spectral variations are modelled in terms of direct light from 
the active nucleus and time-delayed reprocessed light from surrounding gas clouds. We 
use the photoionisation code CLOUDY to evaluate line and continuum reprocessing 
efficiencies e(A, $h, -Wh, &) for clouds of hydrogen density jih and column density 
TVh exposed to hydrogen-ionising photon flux $h- The gas distribution is described by 
a 5-dimensional map, the differential covering fraction f(R, 9, nn, N-r, v), which we re- 
construct from the 2-dimensional data -Fa (A, t) by using maximum entropy techniques. 
Tests with simulated data and a variety of geometries (shells, rings, disks, clouds, 
jets) are presented to illustrate some of the capabilities and limitations of the method. 
Specifically, we reconstruct 3-dimensional geometry-density maps, f(R,9,nn), by fit- 
ting to well-sampled light curves for the continuum and 7 ultraviolet emission lines. 
The maps are distorted in ways that we understand and discuss. The most successful 
test recovers a hollow shell geometry, determining correctly its radius and density. The 
data constrain the ionisation parameter U oc i? -2 ^ 1 to about 0.1 dex, the radius R to 
0.15 dex, and «h to 0.3 dex. We expect better constraints to arise from future fits us- 
ing more lines and velocity profiles as well as velocity-integrated line fluxes. The maps 
are sensitive to the assumed distance, offering some prospects for using emission-line 
reverberations to measure the luminosities and distances of active galactic nuclei. 

Key words: methods: data analysis - techniques: spectroscopic - quasars: emission 
lines - galaxies: Seyfert - galaxies: active - ultraviolet: galaxies 



1 INTRODUCTION 

1.1 Photoionisation Modelling 

Our understanding of the emission-line spectra of active 
galactic nuclei (AGNs) has progressed considerably since the 
late 1970s, when photoionisation calculations were first used 
to model observed spectra in an effort to determine the phys- 
ical characteristics of the emission-line gas (e.g., Davidson & 
Netzer 1979). Sophisticated photoionisation codes now em- 
body our current detailed knowledge and understanding of 
atomic physics and physical processes in astrophysical plas- 
mas. The evolution of these codes toward increasing realism 



is guided by detailed comparison of observed and predicted 
spectra, made possible by advances in high-speed comput- 
ing. 

In the photoionisation code CLOUDY (Ferland et al. 
1998), a prescribed spectrum L\ is directed into a gas cloud 
at radius R characterized by its hydrogen density nn and 
hydrogen column density Nu- The temperature and ioni- 
sation state of the gas in successive zones inside the cloud 
are evaluated, as well as the line and continuum fluxes that 
emerge from the front and back faces of the cloud. For given 
elemental abundances, the calculated emission-line flux ra- 
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tios depend upon L\, R, nu, and Nu- However, the primary 
dependence is upon the ionisation parameter 

Qh = _£h_ 

4irR 2 nuc tihc 
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is the rate at which the source emits hydrogen-ionising pho- 
tons (A < A H = 912A), and 



$ H = 



Qh 



te& (3) 
is the flux of hydrogen-ionising photons incident on the 
cloud. By comparing observed emission-line ratios (which 
are rather similar for most objects) with predictions from 
single-cloud photoionisation calculations, "typical" physi- 
cal conditions in the emission-line gas are inferred to be 
nu ~ 10 9 ' 5 cm~ 3 , U ~ 1(T 2 , $ H ~ 10 18 cm- 2 s _1 (Davidson 
& Netzer 1979). For these parameters the distance of the 
gas cloud from the nucleus is 
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where L44 ~ Qh hc/Xu is the hydrogen-ionising luminosity 
of the source in units of 10 44 erg s _1 . 

Since the derived gas temperatures are low (T ~ 10 4 K), 
the broad emission-line profiles (Vfwhm ~ 5000 km s _1 ) are 
interpreted as arising from bulk motion of the emission-line 
gas, most likely arising from gravitational acceleration of the 
gas by the nucleus. A rough virial mass may then be esti- 
mated by combining the Vfwhm with the radius estimated 
above, 



M~4x 10 S M, 



Vfwhm 



5000 km s 



R 



100 light days 



(5) 



The rough similarity of observed emission-line ratios in 
high- luminosity quasars and lower-luminosity Seyfert galax- 
ies suggests that although these active galactic nuclei span 
a wide range of luminosity, their emission-line regions are 
nevertheless characterized by roughly the same ionisation 
parameter. Moreover, in many cases comparison of veloc- 
ity widths and profiles of different emission lines suggests 
that U is roughly independent of R and hence iih oc R~ 2 . 
These regularities in the observed properties of quasar emis- 
sion lines have become recognized as a problem because it is 
not known what mechanisms may regulate conditions in the 
emission-line gas to ensure a constant ionisation parameter. 

A plausible solution to this "fine-tuning" problem has 
recently emerged. The LOC model (Locally Optimally- 
emitting Clouds, e.g., Baldwin et al. 1995 ) proposes that the 
photoionised gas is not adequately characterized by a single 
density and column density at each radius, nu(R), Nu(R)- 
Instead, a whole population of gas clouds, f(R, nu, Nu), pro- 
vides for a wide range of nu and Nu at each 7?. The total 
emission in a given line is then obtained by summing con- 
tributions from all these different cloud types, but most of 
the emission arises from the specific cloud types that are 
maximally efficient in producing that particular line. The 
photoionisation models indicate that high reprocessing ef- 
ficiency in a given line is achieved by a sub-population of 



the clouds defined primarily by a specific range of ionisa- 
tion parameter. In the LOC model, the optimal conditions 
are obtained at many radii in a given source by employing 
lower density gas when farther from the nucleus, and over a 
wide range of source luminosity by using gas at larger radius 
in higher- luminosity sources. Thus the apparent fine tuning 
is a consequence of photoionisation physics rather than spe- 
cific conditions in the emission-line gas. 



1.2 Reverberation Mapping 

The assumption that the emission lines are driven by pho- 
toionisation is amply validated by variability monitoring 
campaigns. The active nuclei in half a dozen Seyfert 1 galax- 
ies have been studied by the AGN Watch collaboration (Al- 
loin et al. 1994), combining IUE and ground-based obser- 
vations to monitor variability in the ultraviolet and opti- 
cal spectra. The AGN Watch campaigns extend over many 
months, with time sampling typically a few days (Peterson 
1993). In all cases the observed variations in the contin- 
uum flux, F c (t), are accompanied by correlated variations 
in emission-line fluxes, F((t). Moreover, the emission line 
variations lag behind those seen in the continuum. The time 
delay, r, can be interpreted as light-travel time. This effec- 
tively measures the radial distance, 



R. 



(1 



(6) 



from the nucleus (more accurately from the continuum- 
emitting regions close to the nucleus) to the line-emitting 
region. Here z is the cosmological redshift of the source. 

The time delays are usually estimated by cross- 
correlating the observed line and continuum light curves 
(e.g., Edelson & Krolik 1988, Koratkar & Gaskell 1991 
White & Peterson 1994). These results have revealed time 
delays of days to weeks, with shorter delays in the higher 
ionisation lines, consistent with more highly ionised gas oc- 
curring closer to the nucleus (e.g., Clavel et al. 1991). 

More detailed information can be extracted by "rever- 
beration mapping" (Bahcall, Kozlovsky & Salpeter 1972, 
Blandford & McKee 1982). A delay map $(t) can be recov- 
ered by fitting to high-quality densely-sampled light curves, 
for example by using the linear echo model 



#(r) F c (t-r) dr . 



(7) 



Maximum entropy fitting techniques have proven to be use- 
ful in this regard (e.g., Home, Welsh & Peterson 1991, Krolik 
et al. 1991, Home 1994), as have regularized linear inversion 
methods (Vio, Home & Wamsteker 1994, Pijpers & Wanders 
1994, Krolik & Done 1995). 

For the strongest lines, C IV in particular, data qual- 
ity can be sufficient to record reverberation effects in the 
velocity profiles, allowing the recovery of kinematic infor- 
mation in the form of a velocity-delay map ^(v,t) (e.g., 
Wanders et al. 1995, Ulrich & Home 1996, Done & Krolik 
1996). Results to date show that the red and blue wings of 
the Civ emission line respond together with smaller delays 
than the line centre. Models involving purely radial inward 
or outward motions are ruled out, but a variety of models 
involving primarily random or azimuthal motions, as in a 
disk, are permitted. 
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Time delays measured from H/3 and optical contin- 
uum light curves of AGNs with a range of luminosities 
(10 41 ' 8 erg s" 1 < AL A (5100A) < 10 45 ' 7 erg s" 1 ) establish 
an empirical relationship of the form 

(Kaspi et al. 2000). This is consistent with the photoionised 
region increasing in size as the ionising luminosity increases. 

The echo-mapping experiments have forced us to re- 
alize that emission-line gas is present much closer to the 
nucleus than the radius estimated in Equation (4) on the 
basis of single-cloud photoionisation models. For a given U, 
the smaller radius requires a higher density, riH ~ 10 n cm~ 3 
(Ferland et al. 1992). In earlier work, the ClIl]/ClV ratio 
was used to argue against such high densities. This argu- 
ment may be incorrect, however, if Cm] and Civ form in 
different regions with different densities, as is suggested in 
several cases where the C ill] lag is longer than the C IV lag 
(Clavel et al. 1991, Reichert et al. 1994). The smaller radii 
also reduce the inferred virial masses. 



1.3 Quasars as Cosmological Probes 

Quasars are potentially important as cosmological probes. 
Their high luminosity renders them visible at large redshifts. 
However, quasars are poor standard candles, with luminosi- 
ties ranging over many orders of magnitude. In low-redshift 
samples the emission-line ratios and equivalent widths are 
correlated with luminosity (Baldwin 1977; Baldwin et al. 
1978; Kinney et al. 1990), but the dispersion is too large 
to be of much use in estimating luminosities or distances. 
Separating luminosity and evolution effects is also difficult 
with magnitude-limited samples that tend to include high- 
luminosity objects at high redshift and low-luminosity ob- 
jects at low redshift. 

Distances may alternatively be estimated by using light 
travel time delays from reverberation effects to measure the 
linear size of something whose angular size can be directly 
observed or inferred from the observed spectrum. For exam- 
ple, the distance to supernova 1987A was estimated from 
time-delayed enhancements in photoionised emission lines 
produced when the ultraviolet flash from the explosion first 
reached the inclined circumstellar ring that is resolved by 
HST (Panagia et al. 1991). Can we apply a similar method 
to AGNs? 

For AGNs a method based on reverberations in con- 
tinuum emission from the irradiated surface of the accre- 
tion disc has been demonstrated by Collier et al. (1999). 
A steady-state disc is predicted to have a characteristic 
temperature profile of the form T oc R~ 3 ^ 4 . With r ~ 
(1 + z)R/c from light travel time, and kT ~ hc(l + z)/4\ 
from blackbody radiation, the temperature decreasing with 
radius implies a time delay increasing with wavelength as 
t oc A 4 / 3 (l + z)^ 1 ^ 3 . Just such an effect was found, with 
delays of order 1-2 days between ultraviolet and optical con- 
tinuum variations in the Seyfert 1 galaxy NGC 7 469. The r e- 
sulting redshift-independent distance yields Ho \J cos i/0.7 = 
42 ± 8 km s _1 Mpc _1 , where i is the inclination of the disk 
axis to the line of sight, expected to be less than 60° for 
Seyfert 1 galaxies. 



In §3 of this paper we discuss a new method of esti- 
mating distances to AGNs based on reverberation in the 
photoionised emission lines. The basic idea is to require the 
radii of photoionised regions estimated from time delays to 
match the radii derived from comparison of photoionisation 
calculations with the observed emission-line ratios. This may 
provide a new direct method of determining distances, based 
on a straightforward interpretation of time delays and pho- 
toionisation physics. 



1.4 Quasar Tomography 

The goal of quasar tomography, as developed in this paper, is 
to unify echo mapping and photoionisation modelling. Echo 
mapping reveals the sizes of emission-line regions while pho- 
toionisation modelling uses emission-line flux ratios to deter- 
mine physical conditions within the emission-line gas. Cur- 
rent echo mapping techniques define delay maps for each 
line independently of the other lines and without regard 
to photoionisation physics. More complete information can 
be extracted from high-quality time-resolved emission-line 
spectra by combining these two approaches. 

Following the spirit of the LOC model, we offer the 
emission-line gas freedom to fill out a very general 5- 
dimensional map, f(R,0,nn,NH,v), giving the differential 
covering fraction of the gas clouds. We then expect observa- 
tions of time- dependent spectra -Fa (A, i) to reveal the geome- 
try, physical conditions, and kinematics of the gas by placing 
constraints that define the structure of the 5-D cloud map. 

We are aiming to fit simultaneously, in far greater de- 
tail than has hitherto been attempted, the time-dependent 
fluxes and velocity profiles of numerous emission lines 
recorded in high quality spectra. Our methods and assump- 
tions in modelling the reverberating spectrum F\(\, t) are 
developed in Appendix A. Appendix B then discusses the 
maximum entropy techniques we employ to recover the 
cloud map f(R,6,nn,Nn,v) by fitting to the observations 
of -F\ (A, t). 

In §2 we illustrate some of the capabilities and limita- 
tions of quasar tomography by presenting the results of test 
reconstructions from simulated datasets. A variety of possi- 
ble geometries is considered. In §3 we assess prospects for 
using this approach to derive redshift-independent distances 
and luminosities. Concluding remarks are made in §4. 



2 CLOUD MAPS RECOVERED FROM 
SIMULATED DATASETS 

2.1 Appearance of Different Geometries 

It is important to realize that most types of observational 
data obtained from a distant unresolved object are unaf- 
fected if the object is rotated by an angle <j> around the line 
of sight. A possible exception is polarimetric data, where 
the position angle of linearly polarized light is measured. In 
quasar tomography we employ information from time delays 
and emission-line ratios that are sensitive to R and 8 but not 
to <f>. Thus we should not expect to be able to recover the 
3-dimensional source geometry, f(R, 8, <jj), but rather the 2- 
dimensional geometry, 
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6 Geometries for R= 1 2 i=45 
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Figure 1. Maps showing the appearance of six different geome- 
tries after being rotated around the line of sight (X axis). 

f(R,0) = J f(R,6,4>) d4>, (9) 

that arises by rotating the 3-dimensional geometry around 
the line of sight to the observer. 

To illustrate the effect of this ambiguity, Fig. 1 displays 
the appearance of several simple geometries after projection 
around the line of sight onto the X, Y plane. Here X = 
—RcosO is measured from the nucleus towards the observer, 
whose location is off to the right-hand side of the figure, and 

Y — it! sin # is measured perpendicular to the line of sight. 
The ambiguity in <f> means that we can project any given 
geometry by rotating it around the line of sight. This yields 
a projected 2-dimcnsional geometry that is symmetric to 
reflection through the X axis. A single cloud (upper-right 
panel) projects to a spot at the appropriate radius R and 
azimuth 0, plus its mirror image reflected through the X 
axis. A hollow spherical shell (upper-left panel) projects to 
a complete ring with radius R equal to that of the shell. 

An inclined ring or annulus (middle-left panel) projects 
to an arc of radius R that straddles the +Y axis, plus its 
mirror image on the —Y axis. The arc spans 9 = 90° ± i, 
where i is the inclination angle between the ring axis and 
the line of sight. The figure shows the projection of a ring in- 
clined by i = 45° . If we tilt the ring axis toward the observer, 
the two arcs become shorter, shrinking to a pair of spots at 

Y — ±R as i — > 0. Tilting the ring axis away from the ob- 
server makes the arcs longer until their endpoints meet to 
form a complete ring as i — > 90° . This edge-on ring may still 
be distinguished from the hollow shell, however, because a 



relatively large fraction of the ring's circumference projects 
to the region near the endpoints of the arc. 

An inclined disk (lower-left panel) is assembled from 
a set of co-axial annuli with different radii. The disk an- 
nuli project to arcs of different radii that fill out the sector 
bounded by the smallest and largest arcs. This makes an 
X-shaped pattern with the sectors above and below the X 
filled in and empty sectors to the sides. As with the inclined 
ring, tilting the disk toward (away from) the line of sight 
closes (opens) the filled sector of the X. 

A 1-sided hollow conical jet (middle-right panel) is eas- 
ily recognizable from its projection to a cone emerging from 
the nucleus. The jet inclination i is the angle between the 
cone and the X axis, and the jet opening angle u is the same 
as that of the cone. Of course the mirror image of the jet, 
projected across the X axis, is also present after rotation 
about the line of sight. If we increase the jet opening angle, 
the corresponding cone becomes wider, eventually overlap- 
ping with its mirror image when uo > i. Finally, twin jets 
extending in opposite directions (lower-right panel) project 
to an X-shaped pattern of cones. 

The main point we wish to make here is that these six 
geometries, and many others, have distinct appearences af- 
ter projection by rotation around the line of sight onto the 
X — Y plane that we aim to recover from the data. In §2.3 
we present and discuss maps of these same six geometries 
reconstructed from simulated datasets. 

2.2 Synthetic Light Curves 

Fig. 2 shows a simulated dataset calculated using the meth- 
ods outlined in Appendix A for the geometrically-thin spher- 
ical shell shown in Fig 1. All clouds within the shell are as- 
sumed to have a gas hydrogen density of na = 10 n cm~ 3 , 
and a hydrogen column density of JVh = 10 23 cm -2 . 

The bottom panel shows the driving light curve, F c (t), 
and the responding emission-line light curves F e (t) are 
shown for 7 different ultraviolet emission lines in the 7 pan- 
els above. The data points in these panels were generated by 
first evaluating the true light curve, and then adding a Gaus- 
sian random variate with zero mean and standard deviation 
as shown by the vertical error bars (smaller than the data 
points plotted) to simulate 3% observational errors. The hor- 
izontal dashed lines show the background fluxes fs(A) for 
each line, included in the model to represent constant line 
emission that is not driven by photoionisation. Light curves 
similar to these would be possible to obtain in one year of 
observation by observing once every two days. 

Delay maps ^(r) are shown in Fig. 2 to the left of 
the corresponding light curve for each emission line I. The 
delay maps for this geometrically-thin spherical shell with 
R — 12 light days extend from to 2R/c — 24 light days. 
The prompt response (at r = 0) arises from the near edge 
of the shell, while the response with maximum delay (r = 
2R(1 + z)/c) arises from the far edge of the shell. When 
the line emission is isotropic (e.g., Cm] in Fig. 2) then the 
delay map for the spherical shell is flat-topped, with equal 
response in each interval of r from to 2R(1 + z)/c. For the 
other lines the delay maps have a positive slope (response 
increasing with increasing r) due to their inward anisotropy 
(enhanced emission at large delays arising from the far edge 
of the shell) . 
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Figure 2. Simulated continuum and emission-line light curves 
for a thin spherical shell. The continuum light curve (lower panel) 
convolved with the appropriate delay map (left panels) gives the 
corresponding emission-line light curve (right panels). The delay 
maps account for light travel time delays as well as non-linear and 
anisotropic emission-line responses to variations in the ionising 
radiation. See text §2.2 for details. 

The exact angular pattern of anisotropic line radiation 
depends of course on details of geometry and radiative trans- 
fer in the reprocessing region. These details are beyond the 
scope of current 1-dimensional photoionisation codes. We 
are therefore handling the anisotropic line emission by us- 
ing CLOUDY to evaluate the inward and outward emission 
of the line, and then interpolating linearly in cos# to ob- 
tain the emissivity in direction 9. For further discussion see 
§A4.1. 

The driving light curve (bottom panel of Fig. 2) varies 
by a factor ~ 4 in this example, and the delay maps can 
change appreciably when the ionising flux changes by such a 
large factor between the bright and faint states. To illustrate 
this, the left panels of Fig. 2 give delay maps corresponding 
to both the bright state (solid curves) and the faint state 
(dashed curves) corresponding to the maximum and mini- 
mum of the driving light curve. A third delay map (dotted 
curve) is the difference between the bright and faint state 
delay maps, scaled to have the same peak as the bright state 
map. This presentation makes it easy to spot non-linear line 
responses. The scaled difference map (dotted curve) will dif- 
fer from the bright state map (solid curve) only if there is a 
significant non-linear response in that emission line. It can 
be seen that the responses in all lines are positive and es- 
sentially linear, with the exceptions of Lya, and Mgll. 



Lya shows a marked non-linear response primarily due 
to a change in the line radiation pattern between the faint 
and bright states. In the faint state, Lya has a strong in- 
ward anisotropy because its large optical depth causes the 
bulk of its emission to emerge from the irradiated faces of 
clouds (Ferland et al. 1992; O'Brien et al. 1994). Thus we 
see stronger Lya emission from clouds on the far side of the 
shell than we do from those on the near side. In contrast, 
Lya emission is almost isotropic in the bright state because 
the hydrogen-ionisation front pushes through the clouds, so 
that those on the near side of the shell become visible. Thus 
in the top left panel of Fig. 2 we see that in the low state 
(dashed delay map) the prompt Lya response at r = 0, 
arising from the near side of the shell, is essentially zero. In 
the high state (solid delay map), the Lya response is almost 
independent of r, indicating nearly isotropic Lya emission 
from the spherical shell. 

The Mgll response is also non-linear. Moreover, the 
Mgll response on the far side of the shell is negative, mean- 
ing that Mgll emission decreases with increasing ionising 
flux. In the faint state, the Mgll radiation is dominated by 
emission from the irradiated sides of the clouds, producing 
a strong inward anisotropy with very low prompt response. 
In the bright state, the Mgll emitting zone becomes sig- 
nificantly depleted as the hydrogen-ionisation front moves 
deeper into the cloud, and the line emerges more isotropi- 
cally. The depletion of the partially ionised zone where this 
line is formed results in a net negative response for this line 
(see also Goad et al. 1993; O'Brien et al. 1995). The nega- 
tive response of Mg II can be seen directly in the light curves, 
where the Mg n flux has minima at times when the other line 
fluxes have maxima. 

In the spherical shell example we are considering, di- 
agnosis of the delay maps reveals isotropic and anisotropic 
emission, linear and non-linear responses, and positive and 
negative responses in the various emission lines. Each of 
these effects conveys information about the possible location 
and physical state of the gas clouds that are responding to 
changes in the flux of ionising radiation represented in the 
dataset by the driving light curve. The information we seek 
is not readily available, however, as it appears in the de- 
lay maps rather than in the observed light curves. We must 
therefore consider to what extent it is possible to recover 
the geometry and physical conditions by fitting to the light 
curves. 

2.3 Fits to Synthetic Light Curves 

Fig. 3 shows our fit to the simulated dataset of Fig. 2. This 
fit is accomplished, using the maximum entropy methods 
discussed in in Appendix B, by adjusting the driving light 
curve L(t), the background spectrum Fb{\), and the cloud 
geometry and density map f(R, 9, hh). We hold the column 
density fixed at the correct value Nn = 10 23 cm~ 2 . Note that 
our model does not assume a spherical shell geometry with 
all clouds having the same density, but rather determines 
what geometry and density distribution are required in order 
to fit the light curves. 

The predicted light curves are required to achieve a 
good fit, as judged by the criterion x 2 /N = 1, where 
N = 968 is the number of data points. The fitted model 
predicts line fluxes that are not quite as high as the data 
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Figure 3. A fit to simulated continuum and emission-line light 
curves for a thin spherical shell. See text §2.3 for details. 



in the peaks of the light curve, particularly for Lya. This 
is the principal defect in the fit. It arises because our maxi- 
mum entropy fit seeks the "smoothest" functions that fit the 
data points. The reconstructed delay maps, \&^(r), are noisy 
at the 10% level, but they do bear a satisfying resemblance 
to the true delay maps shown in Fig. 2. 

The maximum entropy fitting techniques are discussed 
in detail in Appendix B, but it may be appropriate here 
to summarise the method and discuss a few issues. The fit 
arises by iterating from initial guesses for the driving light 
curve L(t), the background spectrum Fb(X), and the cloud 
map f(R,d,nii). The iteration seeks to improve the fit to 
the data, aiming for \ 2 /N = 1, while also maximizing the 
"entropy" to keep the model "as simple as possible", in a 
well defined sense. The iteration either converges or else 
informs you that it has failed to converge. The converged 
model represents the "simplest" model that fits the data 
with x 2 /N = 1. 

The maximum entropy fit is unique in the sense that it 
does not depend on the initial guesses, and can be reached 
from a wide range of starting points. However, in a differ- 
ent sense the maximum entropy fit is not unique because 
there are a variety of ways to define what you mean by "as 
simple as possible". For the fits presented in this section, 
we follow our usual practice of defining "simple" to mean 
"smooth" . Maximizing the entropy "steers" each parameter 
of the model toward values in adjacent pixels. This delivers 
the "smoothest" background spectrum, driving light curve, 
and cloud map that succeed in fitting the data. When two 



log n H = 23.00 log * H = 21.7 : 19,7 , 22.4 : 20.3 




Figure 4. Maps derived from a maximum entropy fit to the sim- 
ulated data shown in Fig. 3. The reconstructed geometry f(R, 8), 
shown in the lower-right panel, reveals the basic form of a hollow 
spherical shell, with some smearing along the iso-dclay parabolas, 
r oc R(l + cos 9), opening to the right. The other panels display 
maps of the responses of different emission lines to increases in 
ionising radiation from the central source. The greyscale assigns 
white and black to the minimum and maximum responses oc- 
curing in each line. Positive responses arc exhibited by all lines 
except Mgll, which has a negative response on the far side and a 
positive response on the near side of the shell. See text §2.3.1 for 
further discussion. 

models fit the data equally well, the smoother model is pre- 
ferred. When the data require a sharp feature in the map, 
e.g. a peak in the spectrum or light curve, or clouds con- 
centrated at some radius or with some density, maximizing 
entropy spreads out the feature as much as possible in any 
directions permitted by the data. In this way the resolution 
of the map is determined by the quality of the data, and 
the sensitivity of the data to each part of the map. We will 
point out examples of this effect below. 

2.3.1 Density-Geometry Maps of the Spherical Shell 

Our maximum entropy fit to the light curves in Fig. 3 re- 
covers a 3-dimensional map, f(R,6,nn), specifying the ge- 
ometry and density distribution of the emission-line clouds. 
We visualize this by displaying 2-dimensional projections, 
f(R,8) in Fig. 4 and f(R,na) in Fig. 5, which we discuss 
below. 

Fig. 4 shows the reconstructed geometry. The lower- 
right panel gives the differential covering fraction, f(R,9), 
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Figure 5. The density-radius map f(R,nu) recovered from the 
fit to the simulated data shown in Fig. 3. The correct density 
«H = 10 11 cm -3 is recovered, with ambiguity along the direction 
of constant ionisation parameter U oc R~ 2 n^ { 1 . 



and the other panels show the responses in the seven ultra- 
violet emission lines. The basic morphology of a thin spher- 
ical shell, with radius R/c ~ 12 light days, is clearly recog- 
nizable, though distorted in ways that we understand and 
discuss below. The inward anisotropy in most of the lines 
(stronger response from the far side of the shell) is recov- 
ered. Mgn looks peculiar in the plot because clouds on the 
far side of the shell have a negative response in this line. A 
hollow inner region, R 10 light days, is clearly recovered. 
On the far side of the shell, the cloud distribution is con- 
fined to within £ 2 light days of the correct radius. At other 
azimuths, however, the clouds have spread outward along 
iso-delay parabolas, which open toward the right. Striations 
along the iso-delay parabolas are visible mainly at the top 
and bottom of the shell. 

Fig. 5 shows the density-radius projection, f(R,nn). 
Here we see that the correct density, n ~ 10 n cm -3 , is recov- 
ered, though with smearing along the direction of constant 
ionisation parameter, U oc R~ 2 n^ . The line ratios evidently 
constrain U to within about 0.1 dex, while the constraints 
on 71h and hence R are weaker, 0.3 and 0.15 dex respec- 
tively. The data allow clouds can move to a larger (smaller) 
R provided they reduce (increase) iih to maintain a fixed 
ionisation parameter, i.e. tin oc R~ 2 . The line ratios con- 
strain iih because each emission line becomes optically thin, 
due to thermalisation or de-excitation, above a different crit- 
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Figure 6. Maps of the six test geometries illustrated in Fig. 1 
reconstructed from simulated light curves with time sampling and 
signal-to-noise levels similar to those shown in Fig. 3. 

ical density (e.g. Hamann et al. 2002, see also Fig. 5). This 
inhibits inward more than outward changes in R. 

Motion of clouds along iso-delay parabolas preserves the 
time-delay constraints imposed by the light curves. On the 
near side, clouds constrained by r can move in R with only 
a little change in 9, while on the far side they can move in 
6 more than R. This explains why the reconstructed shell 
in Fig. 4 is sharper on the far side than on the near side. 
Motion in 6 must also be inhibited to some extent by the 
anisotropic radiation patterns in different lines. 



2.3.2 Smooth Maps of Six Geometries 

Encouraged by the results obtained for the spherical shell 
geometry, we decided to more thoroughly test the capabili- 
ties of quasar tomography by reconstructing maps from sim- 
ulated data for a variety of test geometries. Fig. 6 exhibits 
our reconstructed maps for the six test geometries illustrated 
in Fig. 1. The maps in Fig. 6 represent the "smoothest" pos- 
itive maps that fit the data with x 2 /N = 1. 

Comparison of Figs. 1 and 6 indicates that significant 
ambiguities in the geometry remain after imposing con- 
straints from the ultraviolet emission-line light curves. In 
all six cases the distorted geometry arises because the data 
constrain U and r more tightly than iih, R or 9, as we dis- 
cussed above in some detail for the hollow shell geometry. 
While the hollow shell geometry is recognizable in the map, 
the map of the inclined ring looks quite similar. The map 
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Figure 7. As in Fig. 6 except that here the maps are recon- 
structed with a preference for front-back symmetry that is char- 
acteristic of point-symmetric and axi-symmetric geometries. 



of the single cloud is greatly extended along the iso-delay 
parabola, making it appear similar to the map of the 1-sided 
jet. The 2-sided jet is hardly recognizable, and the disk is to- 
tally unrecognizable, presumably because these geometries 
extend over large ranges in R, 6, and r. 

The present tests show that the reconstructed geom- 
etry can be significantly distorted while still fitting the 
lightcurves. Nevertheless, we consider that even distorted 
maps such as these would represent interesting progress 
given our current limited understanding of the inner regions 
of active galactic nuclei. Furthermore, ambiguities should 
decrease as more lines are included in the analysis, since 
each line responds with a different sensitivity to U and nu. 
We also expect to reduce the ambiguity by fitting velocity- 
resolved rather than velocity-integrated lightcurves, since 
the line ratio constraints are then available separately for 
the clouds in each velocity bin, rather than only for the 
sum of all the clouds. We intend to test these conjectures 
in future work. We consider next how to make use of prior 
information on the symmetry of the geometry. 



2.3.3 Maps Steered toward Front-Back Symmetry 

The maximum entropy mapping formalism (MEM), devel- 
oped in Appendix B, allows us to incorporate prior infor- 
mation by "steering" the fit toward models that obey some 
type of symmetry. This is done by defining the entropy in 
such a way that a local maximum occurs when the model is 



symmetric. MEM then selects the model that fits the data 
with x 2 /N = 1, and is also "simple" in the sense of being 
"as close as possible" to the desired symmetry. 

At present we have relatively little prior information on 
the distributions and properties of the emission-line clouds 
in AGNs, though this may change in the future. If we had 
good reason to assume that the distributions were spher- 
ically symmetric, for example, we would be able to make 
better maps by making use of rather than ignoring that prior 
information. The maps presented so far are steered toward 
blurred copies of themselves, so that the only preferred sym- 
metry is a local smoothness. We now consider how to give 
preference to geometries with specific symmetries. 

It may be helpful to think about this "steering" as a 
way of asking different questions in order to explore a range 
of maps that are consistent with the data. For example, you 
may wish to ask "Are the data consistent with a spherical 
geometry?" MEM delivers an answer to this question in the 
following way. If there is a spherically-symmetric map that 
fits the data with \ 2 / N = 1 , MEM finds that map. Usually 
there will be many spherical geometries, with different radial 
profiles, that fit the data. In that case MEM finds the one 
that is "as smooth as possible" in the radial direction. If 
no spherical geometries are consistent with the data, MEM 
finds the map that achieves \ 2 /N = 1 and is "as close as 
possible" to spherical symmetry. 

Steering toward spherical symmetry is very appropri- 
ate for the hollow shell geometry, and could considerably 
reduce the distortions we noted in our reconstructed map. 
However, spherical symmetry is not so appropriate for the 
other five geometries illustrated in Fig. 1. We therefore opt 
here for a less restrictive symmetry. As discussed in §2.1, 
any point-symmetric or axi-symmetric geometry has a map 
that is symmetric in X, since for each cloud on the near side 
there is a cloud in the corresponding zone on the far side of 
the nucleus. We note that four of the six geometries in Fig. 1 
are symmetric to rotation about an axis that is inclined to 
the line of sight. When these axi-symmetric geometries are 
rotated around the line of sight, the result is symmetric be- 
tween the near and far side, f(X,Y) — f(—X,Y). We may 
therefore expect to be able to "improve" our reconstruc- 
tions of these axi-symmetric geometries by "steering" them 
toward front-back symmetry. The other two maps, which 
don't obey this symmetry, will be altered by this inappro- 
priate steering, but perhaps not by enough to matter. 

Fig 7 exhibits the maps that arise when preference is 
given to maps with front-back symmetry. All six test ge- 
ometries are now rather more easily recognizable in the re- 
covered maps. Smearing is still evident along the parabolic 
iso-delay surfaces, but this effect is reduced by the front- 
back symmetrizing in comparison with the results in Fig 6. 
Be reminded that in all cases the maps fit the data with 
\ 2 /N = 1. Differences among the corresponding maps in 
Figs. 1, 6 and 7 illustrate ambiguity in the geometry that is 
permitted by the data. The basic resolution is set by the time 
delay accuracy, but ambiguity arises primarily because the 
density is not well pinned down, allowing the cloud distribu- 
tion to spread in radius along iso-delay parabolas while ad- 
justing the density to hold the ionisation parameter roughly 
constant. We suspect that this ambiguity can be reduced if 
the dataset includes additional line fluxes that place tighter 
constraints on gas densities. 
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We now briefly discuss each of the maps with a view 
toward understanding what might be learned about the ge- 
ometry by interpreting such maps. The best reconstruction 
is that of the hollow shell, which is quite clearly defined 
by its map in Fig 7. The tight radial resolution (~ 2 light 
days) that was achieved only on the back side of the ring 
(— X) in Fig 6, is now transferred to the front side (+X) 
as well. Outward smearing along time-delay parabolas is re- 
duced and modified by the front-back symmetrizing, leaving 
the largest radial smearing (~ 5 light days) near the Y axis. 
Extrapolating these results, as discussed above, we expect 
that this hollow shell geometry could be very accurately re- 
constructed if we steered the model toward spherical sym- 
metry. 

The inclined ring should ideally appear as two arcs 
spanning the Y axis (see Fig. 1). Its map in Fig. 7 is su- 
perficially rather similar to that of the hollow shell. The 
two geometries could probably be distinguished, however, 
based on the azimuthal distribution of the gas around the 
ring. In particular, the inclined ring is deficient in clouds 
near the X axis, and has four bright points at (X, Y) = 
(±i?sini, ±_Rcosi), corresponding to the nearest and far- 
thest points on the inclined ring. These features are clearly 
visible in the map and would permit an estimate of the ring 
inclination. Like the hollow shell map, outward radial smear- 
ing is expected and evident where the arcs cross the Y axis. 

The inclined disk and twin jet maps in Fig. 7 are closer 
to their correct geometries, shown in Fig. 1, but are still 
not very easily recognisable. Both maps do indicate cloud 
distributions spanning a wide range of radius, in striking 
contrast to the confined radial range found in the hollow 
shell and inclined ring maps. The inclined disk map has less 
gas along the X axis and more along the Y axis, but the 
disk geometry and inclination are not clearly identifiable. If 
you were to assume a disk geometry, the map implies an 
inclination closer to 45° than to 0° (all clouds on Y axis) or 
90° (all clouds on X axis). 

The twin-jet map has a very clear deficit of gas in a wide 
zone around the Y axis, with a fairly well defined edge. How- 
ever, the two jets blend together across the X axis. Some- 
what better resolution would be required to infer the correct 
geometry. 

The single cloud and 1-sided jet geometries have clouds 
exclusively on the near side of the nucleus, with no corre- 
sponding clouds on the far side. This violates the preferred 
front-back symmetry employed in reconstructing the maps 
in Fig 7. Maximizing the entropy with the preferred symme- 
try encourages these maps to develop spurious far-side copies 
of their near-side features. However, it seems that even in 
these highly asymmetric cases the maps are not much af- 
fected, and are quite similar to or perhaps somewhat better 
than those in Fig 6. The quality of the light curve data 
constraints must be sufficient in this example to rule out re- 
sponses at larger time delays, thus over-riding the preferred 
symmetry. 



3 QUASAR DISTANCES 

So far we have assumed that the distance of the active galaxy 
is known. If the distance is incorrect, the assumed luminosity 
will be incorrect and the resulting maps will be distorted. 



How sensitive are the maps to the adopted distance? To 
put this question the other way, can we use our maps to 
determine distances and luminosities of active galaxies? If 
so, then emission-line reverberation studies could provide a 
new method to measure cosmological parameters such as Ho 
and go- 

The observed flux from an active nucleus at redshift z 

is 

L(A e ) _ he Q H S(A e )C(i) 



F(X) = XFx(X) 



Ah 



where L(\) = XL\(X) is the luminosity, A, 
the emitted wavelength, 



cz 
Wo 



1 + - 



(i - qo)z 



(10) 
= A/(l + z) is 

(11) 



1 + qoz + ^1 + Iqoz _ 

is the luminosity distance, where Ho is the Hubble constant 
and qo is the deceleration parameter, Qu is the emission rate 
of hydrogen-ionising photons, 

A H L(A) AhL(A) 



S(X) = 



hc Qli Io XHL W dX 



(12) 



is the dimensionless shape of the photon spectrum, and £(i) 
is a dimensionless continuum anisotropy factor, allowing the 
spectrum to depend on the observer's viewing angle i. 

As we have seen, by using photoionisation models to fit 
measured emission-line flux ratios we can in principle con- 
strain the hydrogen density nu and the ionisation parameter 
U = $h/(«hc), where $h is the hydrogen-ionising flux inci- 
dent on the emission-line gas. Since $h = Qh/4ttR 2 , we can 
estimate the "photoionisation radius" 

^.f^r.aC^^^. (13) 



V47T<I>H / 



V fec$ H S(A)C(i) 



Echo mapping experiments give an independent measure- 
ment of the radius, "reverberation radius" , 

C T 



R T 



(14) 



(l + z)(l + cos0) ' 

which is derived from time delay measurements. Equating 
Rq and R T then gives the luminosity distance as 



hc$HS(X e )C(i) 



1/2 



AhF(A) J (l + z)(l+cos< 
The Hubble constant (for small z) is then 



Ho 



cz 



XhF(X) 



hc$HS(x e )<;(i) 



1/2 



z(l + COS0) 



(15) 



(16) 



3.1 Distance— Geometry Correlations 

Most of the quantities on the right-hand side of the above 
expressions for Dl and Ho are constrained by observations: 
XF\ and z are directly observed, <&n is determined from 
the emission-line flux ratios, r is determined from the light 
curve time delays. The spectral shape S(X) may be esti- 
mated from the spectral energy distribution obtained from 
multi- wavelength observations extending from X-rays to the 
infrared. The extreme ultraviolet range is not directly ob- 
servable, and this introduces some uncertainty. Extinction 
and reddening by dust in the host galaxy and in our own 
Milky Way must also be taken into account, of course. 
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The possibility of anisotropic ionising radiation, 7^ 
l,may be more difficult to diagnose. Unification scenar- 
ios suggest that we view Seyfert f galaxies near the pole, 
i 60° , and Seyfert 2 galaxies at higher inclinations so that 
a thick dusty torus obscures our view of the broad emission 
line region. Large equivalent widths of emission lines would 
imply that the photoionised clouds see a stronger continuum 
than we do. This effect may be useful in estimating £, which 
appears to be similar in different objects, judging from the 
similarity of their emission-line equivalent widths. 

Ambiguity arises also from cos 9. We expect ( cos 9 ) = 
for any point- or axi-symmetric geometry. However, gas 
could be located preferentially on the far side or near side 
of the nucleus. Very good light curves may constrain the 
shape of the time delay map 'I'(t), rather than just the mean 
value of r, and this would constrain cos 9. The value of cos 9 
is also constrained by the different anisotropies exhibited by 
different emission lines. Thus there would appear to be some 
prospect for constraining cos 9, although we do not expect 
this constraint to be very strong unless the data are very 
good. 

Considering the uncertainty in cost 3 , we expect an am- 
biguity to arise between Ho and the inferred geometry. A 
large value of Ho (or equivalently a lower source luminos- 
ity) may be accommodated by decreasing R and increasing 
cos 9, i.e. moving gas inward and placing it on the far side of 
the nucleus. Because f + cos 9 can be at most 2, this degree 
of freedom should allow up to a factor of 2 increase in the 
inferred value of Ho. In practice this can only be realised for 
an extreme geometry with all gas located on the far side of 
the nucleus. Conversely, a small value of Ho can be accom- 
modated by increasing R and decreasing cos 9, thus moving 
gas toward the observer and out to larger radii so that once 
again, r is roughly preserved. 



Hollow Shell mops for i/ =40,60,90 km s 1 Mpc' 
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Figure 8. Maps of a spherical shell with R/c = 12d and n = 
10 11 cm — 3 recovered from fits to the simulated data shown in 
Fig. 3. The Hubble constant used to generate the fake data was 
Ho = 60 km s -1 Mpc -1 , while the three reconstructions assume 
Hq = 40, 60, 90 km s^Mpc" 1 . 



3.2 H a from Simulated Data 

To investigate prospects for determining Ho, and to ver- 
ify the expected correlations outlined above, we performed 
tests with simulated light curves. We generated fake light 
curves assuming Ho = 60 km s -1 Mpc -1 , and then recon- 
structed the gas distribution using various values of Ho to 
see how that affected the fit to the data and the resulting 
cloud maps. When we assumed a value of Ho fairly close to 
the correct value, it was possible to achieve a good fit to the 
light curve, as judged by /N — 1, but with a distorted 
geometry. When we assumed a very discrepant value of Ho, 
it became impossible to achieve X*/N = 1. By using the en- 
tropy to quantify the distorted geometry, we could identify 
the correct value of Ho as that giving the best fit with the 
simplest geometry. 

Representative results are presented in Fig. 8 for a 
geometrically thin spherical shell geometry, with R/c = 
12 light days, hydrogen density 71h = 10 11 cm~ 3 , and hy- 
drogen column density Nu = 10 23 cm -2 . The geometry 
and density maps shown for Ho = 40, 60, 90 km s _1 Mpc -1 
are reconstructed by fitting to light curves generated for 
Ho = 60 km s _1 Mpc -1 . In all three cases the fit achieves 
X 2 /N — 1. The thin-shell geometry is fairly accurately re- 
covered when we assume Ho = 60 km s _1 Mpc -1 . For 
Hq = 90 km s _1 Mpc -1 , the map is highly distorted with 
most of the gas from the shell displaced to smaller R and 



away from the observer to form a clump on the far side of the 
nucleus. For Ho = 40 km s -1 Mpc -1 , the map is distorted in 
the opposite sense, with gas displaced away from the nucleus 
and toward the observer along iso-delay parabolas. 

The correct value, Ho = 60 km s _1 Mpc -1 , may be rec- 
ognized as the value giving the most symmetric map. The 
uncertainty in Ho resulting from such fits may be quanti- 
fied in terms of the posterior probability distribution derived 
from the fit 

P{H \D) oc P{D\H ) P{H ). (17) 

Here P(Ho) is the prior probability assigned to different 
values of Ho before we consider the data D. Regardless of 
our prior view on the value of Ho, the evidence provided by 
the data D multiplies our prior probability by the factor 

P(D\H ) ocexp{aS*-x 2 /2} • (18) 

This applies for each value of Ho the appropriate penalties 
for failing to fit the data (large x 2 ) an d for deploying exces- 
sive numbers of parameters (large negative aS, where S < 
is the entropy of the map and a > is the mixing parame- 
ter). For more detailed discussion, see Appendix B. 
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4 SUMMARY AND CONCLUDING 
REMARKS 

We have developed a method of reconstructing 5- 
dimensional maps of photoionised emission-line regions by 
fitting predictions of photoionisation models to high-quality 
observations of reverberating emission-line spectra of active 
galactic nuclei. Reverberation effects in the emission lines 
offer information about the (R, 6) geometry and the angu- 
lar emission pattern of the gas clouds, coded as time-delay 
information in the light curves of different lines with re- 
spect to the continuum. At the same time, line profiles give 
the distribution of clouds with line-of-sight velocity v, line 
ratios constrain the density tih and column density Nn of 
the gas clouds, and equivalent widths constrain the covering 
fraction /. Our method brings all these constraints together 
in a global fit to reverberating emission-line spectra in or- 
der to recover a 5-dimensional cloud map, the differential 
covering factor f(R,9,nn,Nn,v). This method may allow 
high-quality observations to reveal the geometry, physical 
conditions, and kinematics of the emission-line gas. 

We use the maximum entropy method to find the "sim- 
plest" positive maps (maximum entropy) that fit the data 
(x 2 /N = !)• The entropy is measured with respect to a 
default map, giving us some freedom to choose what we 
consider to be a "simple" map. We use the entropy to give 
preference to smooth maps rather than noisy ones, and to 
"steer" the map toward various possible symmetries - spher- 
ical, axial, point, etc. The range of maps found for different 
choices of symmetry helps to gauge the extent and nature of 
ambiguities that remain after fitting the data. The Bayesian 
formulation provides posterior probability distributions, and 
Monte Carlo techniques may also be used to assess statis- 
tical uncertainties in the cloud map, the distance, and any 
other parameters that affect the fit. 

Tests with simulated datasets indicate good prospects 
for mapping geometry and density structure from observa- 
tions of emission-line light curves that show evidence for 
time delays. Good results require simultaneous fits to nu- 
merous lines in order to constrain both the ionisation pa- 
rameter and the density. The main ambiguity arises because 
the emission-line lightcurves place tighter constraints on the 
ionisation prameter U and time delay r than on the geome- 
try (R, 9) and physical conditions (JVh,tih). Our simulation 
tests suggest that in some cases the correct geometry is rec- 
ognizable in the reconstructed maps when we give preference 
to front-back symmetry, a property of both point-symmetric 
and axi-symmetric cloud distributions. If geometries are in- 
deed approximately symmetric, then distances can be de- 
rived in order to constrain cosmological parameters, (Ho and 

go). 

The test cases presented in this paper consider only 
a single column density, Nn = 10 23 cm~ 2 . A range of Nn 
at each radius will probably be required to fit observed 
lightcurves (Shields, Ferland, Peterson 1995). Extending the 
photoionisation grid to include a range of Nn is straight- 
forward, but will require more computer time for the iter- 
ative fitting. The additional freedom should also increase 
the ambiguity, but perhaps not too severely because in 
general the emission line spectra of clouds in the range 
10 22 < Ah < 10 24 are less sensitive to Nn than to U and nn 
(e.g. Korista et al. 1997; Goad & Koratkar 1998). 



We have largely omitted to discuss mapping the kine- 
matics of the emission-line gas, focussing instead on map- 
ping the geometry by fitting to velocity-integrated emission- 
line light curves. Extracting kinematic information will in- 
volve simultaneous fitting of light curves at many wave- 
lengths to resolve time variations and hence time delays 
at each velocity in the emission-line profile. Extending our 
analysis to include multiple pixels along the u-axis of the 
cloud map is straightforward. Note that continuum fitting 
and de-blending of lines in the observed spectra will not be 
necessary because the model adds together the predicted 
continuum and the velocity profiles of all the lines. It will 
be interesting to see the extent to which v — R projections 
of f(R,9,nn,Nn,v) provide evidence for virial motions. We 
expect the use of velocity information to reduce the ambigu- 
ities, since constraints will then be available separately for 
many subsets of the cloud population rather than only for 
the sum over all clouds. 

We rely on CLOUDY both to generate our fake datasets 
and to reconstruct from them the cloud maps. In the analy- 
sis of real observations, uncertainties in the atomic data, in 
the elemental abundances in the shape of the ionising con- 
tinuum, and in the assumptions and approximations used in 
the photoionisation calculations will each give rise to errors 
and distortions of the maps. It will therefore be important 
to continue to improve the atomic data and photoionisation 
physics in tandem with attempts to fit CLOUDY predictions 
to observations. 

Finally, we note that the methods developed here for 
interpreting spectral variations in active galactic nuclei can 
also be extended to other types of astronomical objects 
in which an unknown gas distribution is excited by pho- 
toionisation from a compact source of ionising radiation. 
In cases where time-dependent reverberation effects are 
not observed, the observed spectrum alone or, better still, 
spatially-resolved spectral information from combinations 
of narrow-band imaging and long-slit or integral-field spec- 
troscopy can be used as constraints. Maximum entropy fit- 
ting then recovers the "simplest" viable maps of the geome- 
try, kinematics, and physical conditions in the photoionised 
gas. A likely future application is simultaneous fitting of 
narrow-band images and long-slit spectra of nova shells and 
planetary nebulae. 
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APPENDIX A: THE FORWARD PROBLEM 
Al Geometry 

Two parameters, At and t M ax, set the pixel size and radius 
of the region that we wish to map. These need to be cho- 
sen appropriately for each dataset. We recommend setting 
At to the median time step between successive data points, 
and Tmax to a value larger than the largest delay for which 



the data provide evidence of significant response, usually not 
more than a third of the full duration of the observed light 
curves. The model includes a constant background spec- 
trum, §A3.1, to allow for responses on larger timescales. Our 
maximum entropy fits (Appendix B) degrade the resolution 
of the map as much as possible within the constraints set 
by the light curve data. If the pixels are too large, the res- 
olution may be limited by the pixel size rather than by the 
data. If the pixels are too small, the resolution is defined by 
the data, but more computer time is required to accomplish 
the fit. Thus high or low signal-to-noise datasets may require 
smaller or larger At, respectively. 

A 1.1 Shells 

In the model developed for this paper, we divide the spatial 
volume surrounding the nucleus into Nr concentric spherical 
shells with radii 



R(i) = iAR 



(Al) 



where AR = cAt/(l + z) and At is a suitable time interval 
comparable to the time spacing of the observations. The 
time delay range covered by the shells is thus r m i n = to 
Tmax = 2NnAt. Any flux arising from outside this region is 
included in the model in the form of a constant background 
spectrum (§A3), thus neglecting any reverberation effects. 
Note that other partitions are possible, for example equal 
spacing in log R may be more appropriate if a large range of 
radius is to be mapped. 

The exact choice of At is not too critical because in 
fitting the data we first construct predicted light curves uni- 
formly sampled with a time interval At, and then interpolate 
to the actual times of observations. Given an equally spaced 
time series, At can be the time interval between succes- 
sive data points. For unequally-sampled time series, At can 
be the minimum or median time spacing, or some point in 
between. A smaller At may permit resolution of finer struc- 
ture if the data record significant flux changes from one data 
point to the next. A larger At may be adequate if the data 
record no significant flux changes between successive data 
points. 

A 1.2 Zones 



The time delay accounting for light travel time is 
t(R,8) = -(l + z)(l+cos6>) , 



(A2) 



where the angle 6 specifies the direction of the gas cloud 
as seen from the nucleus measured from the direction away 
from the observer. Note that the time delay depends on 6 as 
well as R, as do the reprocessing efficiencies discussed in §A4 
below. We therefore sub-divide shell i into Ng(i) = 2i + 1 
equal- area zones delimited by equal intervals of cos 9. The 
solid angle covered by each zone as viewed from the origin 
is 

47T 

(A3) 



a ^ 47T 



2i + l 

The cosine of our viewing angle for clouds in zone j of shell 
i is then 



cos6< 



\Ne-lJ 



3 



i-1 



(A4) 
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Finally, the time delay is 



t = 2iAt 



\N e - l) 



At(j-l) 



(A5) 



Note that because we choose Ng = 2i + 1, the delay spacing 
is At in all shells. The full range of delays is < r < t max = 
2N R At. 

Because none of the observable data we consider * de- 
pend upon the angle (j> measured around the line of sight, 
we can omit further sub-division of the zones into sectors. 

We wish to emphasize that our use of a shell and zone 
partition of the 3-dimensional volume is completely general 
and does not exclude any class of models. While we em- 
ploy nested shells, this does not imply a restriction to a 
spherically-symmetric geometry, since structure in the an- 
gle 9 is resolvable by the zones that partition each shell. We 
do not partition the zones in the angle <j>, but this does not 
restrict us to models that are symmetric to rotation around 
the line of sight. It is just that we are unable to detect any 
(j> structure that may exist because there are no aspects of 
the data that depend upon the angle (j>- What we can detect, 
then, is the true 3-dimensional geometry f(R, 9, <f>) projected 
by rotation around the line of sight to give a 2-dimensional 
map 



f(R,o) 



I 



f(R,e, 



(A6) 



A 1.3 Cloud covering fraction 

Zone 9 of shell R subtends solid angle AQ and con- 
tains a population of gas clouds that intercepts a fraction 
f(R,9) AR Afl/An of the ionising luminosity. The cloud 
cover arises from different types of clouds characterized by 
their hydrogen density nu, column density Nu, and line-of- 
sight velocity component (Doppler shift) v. The cloud ge- 
ometry is therefore obtained by integrating over the cloud 
types : 



f(R,9) = / f(R,e,nu,N u ,v) dnu dNu dv 



(A7) 



This introduces the 5-dimensional differential covering frac- 
tion map that we aim to reconstruct by fitting to the ob- 
served spectral variations. Since / is a distribution, we can 
use any convenient partition to divide the (R, 9, nu, Nu, v) 
domain. In this paper we employ a partition with equal in- 
tervals of R, cos#, logn H , logAf H , and v. 

The radial dimension of a cloud of density nu and col- 
umn Nu is Ir ~ Nu/nu- To fit within its spherical shell, 
the cloud must satisfy Nu/n H < AR. For AR = 10 15 ' 4 cm 
(1 light day), this requires Nu < 10 26 ' 4 cm~ 2 (nH/10 n cm~ 3 ). 
In this way we may justify excluding clouds from a corner 
of the n H — N K plane, although our present implementation 
does not impose this constraint. 

If y > is a cloud's tangential to radial aspect ratio, 
we may consider clouds that are spherical (y = 1), radially- 
elongated "cigars" (y < 1), or radially- flattened "pancakes" 
(y > 1). The tangential area covered by such a cloud is 



Polarimetry data might provide a means of resolving structure 



A ~ y 2 (Nu/nu) 2 , and the differential covering fraction for 
the ensemble of such clouds is then 

'Nu_\ 2 AR A Q 
.nu J 47T 



f(R,9,n H ,Nu,v) 



y 



n c i dy 



(A8) 



where n c \(R,9,nu, Nu,v,y) is the volume number density 
of clouds. Thus our cloud map f(R,9,nu,Nu,v) does not 
assume spherical clouds but rather represents the combined 
sky coverage of whatever cloud shapes are present within the 
volume. Note, however, that in §A4.1 we adopt an anisotropy 
function that may be more appropriate for spherical clouds 
than for cigars or pancakes. 



A1.4 Cloud shadowing 

Our photoionisation code CLOUDY assumes an unob- 
structed path from the origin out to each cloud. If a cloud is 
present at position (R,9,(f>), self-consistency demands that 
clouds at the same 9, (f> but larger R should be unable to 
"see" the centre. ' For our shell-zone geometry, this con- 
straint can be satisfied if 

r R UAX AO A R 

/(*)=/ f(R,9)^^<l. (A9) 



4tt 

At present we do not implement this as a hard constraint 
upon the solution, but rather we use the entropy to "steer" 
the map toward this constraint (see Appendix B) , and check 
the final result. 



A2 Ionising Radiation 

We assume that the nucleus is a point source emitting a 
time-variable spectrum described by 



L(X,t) = XL x (X,t) = L S(X)L(t) 



(A10) 



Here the time variations of the nucleus are specified by the 
dimensionless light curve L(t), the shape of the nuclear spec- 
trum is given by the dimensionless spectrum 



S(X) 



XuL(X,t) 
f AH L(X,t)dX 



(All) 



and the dimensional luminosity scale factor is Lo ~ 
Quhc/Xu, where Qu is the rate at which the nucleus emits 
hydrogen-ionising photons (A < Ah = 912A) in the reference 
state L(t) — 1. This separable form, with S(X) independent 
of t, assumes spectral shape does not change as the luminos- 
ity varies in time. The assumption can obviously be modified 
if required. The specific shape adopted for the calculations 
in this paper is shown in Fig Al. 

The flux of hydrogen-ionising photons incident on a gas 
cloud located a distance R from the nucleus is 



* H (t) 



QuL(t — t ) 
4ttR 2 



(A12) 



t The ionising spectrum incident upon the clouds is clearly a 
combination of the "bare" spectrum from the nucleus together 
with transmitted, reflected, and diffuse spectral components aris- 
ing from other clouds. Since fully transparent clouds do not con- 
tribute to the overall covering fraction calcuation, covering frac- 
tions of greater than 1 are possible. 
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D L = 



nuclear spectrum adopted for NGC 5548 




log A (A) 



Figure Al. The nuclear spectrum adopted in photoionisation 
modelling of NGC 5548 and in the simulation tests presented in 
this paper. 



Here t is the time at which we receive the reprocessed pho- 
tons from the gas cloud, and the time delay r accounts for 
the longer light travel time on the path from the nucleus to 
the gas cloud to the observer compared with the direct path 
from the nucleus to the observer. 



A3 Time-Dependent Spectra 

The spectrum that we observe at wavelength A and time 
t is modelled as the sum of three components, a time- 
independent background spectrum, Fb(X), time- variable di- 
rect light from the nucleus, Fo(X,t) and reprocessed light, 
Ffl(A,t), arising from gas clouds near the nucleus that are 
responding to ionising radiation from the nucleus: 

F(X,t) = XF x (X,t) = F B (X) + F D (X,t) + F R (X,t) . (A13) 



A3. 1 Background Spectrum 

The background spectrum Fb(X) is included in the model 
to account for sources of light that are effectively constant 
on the timescale spanned by the observed light curves. Such 
sources include starlight from the host galaxy or nuclear 
starburst, and reprocessed light from gas outside the domain 
of the map that respond to the ionising radiation but on 
timescales much longer than the duration of the light curves. 



A3. 2 Direct Light 
The direct light is 



Fd(X, t) — Fq S(A e ) L(t) , 



(A14) 



where S(X) and L(t) are the dimensionless spectral shape 
and time variations, and Fq = Lq/^-kD\ is the dimensional 
flux observed in the reference state with L(t) = 1. Note that 
with the source at redshift z we must evaluate the spectral 
shape at the emitted wavelength A e = A(l + z), and that 
Dl is the luminosity distance, e.g. 



H 



1 + 



(1 - go)f 
1 + q z + VI + 2goi 



(A15) 



where Ho is the Hubble constant and qo is the deceleration 
parameter. 



A3. 3 Reprocessed Light 

The reprocessed light that we see at time t, arises from gas 
at (R, 9) responding to the ionising flux from the nucleus at 
the earlier time t — T, where 

t= ^(l + z)(l + cos6») . (A16) 

The gas cloud therefore sees an incident hydrogen-ionising 
photon flux 



QnL(t - t) 



4nR 2 

We model the spectrum of the reprocessed light as 



(A17) 



F R (X,t) = F S(A, 



o) f 
Jo 



L(t-r)*(T,\\L(t-T))dT . (A18) 



(A19) 



Here Ao = 1215A is a reference wavelength, discussed in §A4 
below. Note that this model accounts for non-linear repro- 
cessing by allowing the delay map *(r, X\L) to depend on 
L, the dimensionless ionising luminosity of the nucleus. Non- 
linear responses arise because the reprocessing efficiencies of 
the gas clouds change when they are exposed to different lev- 
els of ionising flux. The delay map is obtained by summing 
contributions from all types of clouds, each contributing at 
the appropriate delay r, and weighted by the differential 
covering fraction / and reprocessing efficiency e: 

*(r, A|L) = / !§dR drift dNa dv 
x f(R,e,nn,N u ,v) 
x e(A, $h, Tin, JVh, v, 9) 
x 5{t- f (l + z)(l + cos<9)} . 

Here e(A, $h, jih, Nr, v, 6) is the dimensionless reprocessing 
efficiency, discussed below, and the Dirac 5 function ensures 
that the correct time delay is used at each reprocessing site. 



A4 Reprocessing Efficiencies 

The reprocessing efficiencies, e c for the continuum and ee 
for line t, are defined in terms of dimensionless equivalent 
widths, i.e. the reprocessed flux emerging from the cloud di- 
vided by the incident flux XF\ at the reference wavelength 
Ao = 1215 A. These we evaluate with the photoionisation 
code CLOUDY (Ferland et al. 1998) for a pre-specified grid 
in log$H, lognH, and logA?H, for given abundances and 
shape of ionising spectrum (Korista et al. 1997). Linear in- 
terpolation in the grid is then used to evaluate the efficien- 
cies as required. 

The reprocessed light includes both line and continuum 



e(A, $h, Tin, Nft, v, 6) = (X /X)e c (X,$ft,nft,Nft,8) 
+ J2iXoGx(x e )e e {^ft,nH,N H ,e) . 



(A20) 



For the reprocessed continuum light, e c (A, $h, nn, Nft, 9) 
gives the ratio of XF\ for the reprocessed continuum flux 
divided by that for the incident flux at Ao . For the emission 
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Figure A2. The response as a function of density and radius is 
shown for each of the major ultraviolet emission lines for a filled 
spherical volume that reprocesses equal fractions of the ionising 
spectrum in each unit of radius. 



lines ei(^u,nu, Nu,0) gives the flux in line £ emerging in 
direction 9 divided by the incident flux Ao-Fa(Ao)- 

When spectra are required, the line emission is assigned 
a Gaussian distribution in wavelength 



Gx(x t ) 



expj-izf} 



2ttAA 



where 



xe = 



(l + z)X-(l+v/c)X e 
AX 



(A21) 



(A22) 



allows for the Doppler shift by velocity v from the rest wave- 
length Xe- The line width A A represents the spectral reso- 
lution of the data and intrinsic broadening, e.g., Doppler 
broadening due to thermal velocities and the range of veloc- 
ities A« in the velocity bin. 

Fig. A2 examines the responses of the major ultraviolet 
emission lines to changes in the flux of ionising photons. The 
figure shows for each radius and density the difference in the 
line flux between bright and faint states that correspond to 
an increase by 0.8 in log 3>h- In this LOC model, we consider 
a single column density log Ah = 23 cm -2 , and a range of 
densities 7 < logrtH < 14 cm -3 . The radial range considered 
is 2 < R < 24 light days, corresponding roughly to a range 
in the ionising photon flux of 22 > log $h > 20. 

Each emission line reprocesses efficiently in a somewhat 
different region of of the nn-R plane. The regions are delim- 
ited to first order by an appropriate range of the ionisation 



parameter U oc $h/«h oc Qa/(R 2 nn)- For constant tin, the 
reprocessing efficiency diminishes at large $h (small R) be- 
cause the clouds become fully ionised and heat to Compton 
temperatures, and at small $h (large R) because only a 
surface layer on the inward face of the cloud is ionised. Effi- 
cient reprocessing requires a higher density at small R, i.e. 
n u oc Q U /(R 2 U). 

In the LOC model, a range of iin is present at each 7?, 
and so the reprocessing regions can simply move outward 
as the ionising photon flux increases. Clouds on the inner 
edge of the reprocessing region become less efficient while 
those on the outer edge become more efficient reprocessors. 
Lining the inner edge of the reprocessing region is a region of 
negative response, where increasing the ionising photon flux 
decreases the reprocessed line flux (see Goad et al. 1993). 

At constant U, reprocessing efficiencies decrease at large 
riH where collisional processes compete with radiative de- 
excitation of the upper levels of the transitions. This trend 
is seen in all lines except Nv, whose response increases with 
increasing riH- The Nv emission line, for solar abundances 
and the spectral energy distribution assumed here, begins 
thermalising above ~ 10 13 cm -3 . If we extended our grid to 
smaller radius and higher density, the N v emission would 
decrease again, appearing as Civ does above ~ 10 11 cm~ 3 . 



A4-.1 Anisotropic Emission 

The reprocessed radiation is generally emitted with an 
anisotropic angular distribution, typically because emission 
is produced on the inward face of an optically thick cloud 
(Ferland et al. 1992). To allow for anisotropic radiation, the 
inward efficiency ej and outward efficiency eo are computed 
separately. For the angular pattern of the reprocessed light 
we assume 



e(9) = e + (ei - e) cos 9 , 
where the mean efficiency is 
e= (ej + eo)/2 . 



(A23) 



(A24) 



The linear dependence on cos 9 corresponds to isotropic ra- 
diation from each point on the surface of a spherical cloud, 
with the inward and outward fluxes assigned to the inward 
and outward facing hemispheres. 

Fig. A3 shows maps for each of the major ultraviolet 
emission lines for a filled spherical volume that reprocesses 
the same fraction of the ionising radiation in each shell. 
Note that there is a sum over density to obtain the total 
response in each pixel of the map. The sum over density in- 
cludes clouds with both positive and negative response. The 
clouds with positive responses evidently dominate, so that 
in no part of the emission-line region is there a net negative 
response. 

An inward anisotropy (stronger response from the far 
side of the map) is evident for most of the lines. This arises 
from high-density clouds that are optically thick. The ex- 
ception is C ill] , which we treated as isotropic in the present 
calculation. We expect the C ill] inward/total ratio to be 0.5- 
0.6 for clouds that are emissive in this line (O'Brien, Goad 
& Gondhalekar 1994), and we plan to include this small 
anisotropy in future calculations. 
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Figure A3. Maps showing the appearance in each of the major 
ultraviolet emission lines for a filled spherical volume that repro- 
cesses equal fractions of the ionising spectrum in each unit of 
radius. The inward anisotropy in certain lines is evident. 

A5 Delay Maps 

Delay maps $(r) are shown in Fig. A4. Note that each spher- 
ical shell makes a wedge shaped contribution to ^(r) in the 
region < r < 2R/c. The near side of each shell contributes 
at r = while the far side contributes at r = 2R/c. The 
slope of the wedge reflects the degree of anisotropy of the 
line emission, with steeper slopes indicating increased in- 
ward anisotropy. The filled sphere is thus a sum of wedge- 
shaped contributions from the various shells. 

For each emission line the response has been evaluated 
in both a bright and a faint state, and the corresponding de- 
lay maps are shown as solid and dashed curves. Also shown 
as a dotted curve is the difference between the bright and 
faint state delay maps, scaled to the same peak value. The 
dotted curve would be identical to the solid curve if the re- 
sponse is strictly linear, and this is evidently closely satisfied 
for all of the lines. 

The linearity of the responses and the lack of negative 
responses seen in this LOC model are consistent with the 
two implicit assumptions normally made when maximum 
entropy techniques are used to recover delay maps from ob- 
served light curves (e.g., Home, Welsh & Peterson 1991). 

A6 Synthetic Light Curves 

Fig. A5 shows synthetic light curves calculated for the LOC 
model. The bottom panel shows the driving light curve, 
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Figure A4. Delay maps for each of the major ultraviolet emission 
lines for a filled spherical volume that reprocesses equal fractions 
of the ionising spectrum in each unit of radius. Each shell pro- 
duces a wedge of emission extending over < r < 2R/c with a 
slope that increases with the inward anisotropy. 

which was generated assuming that log of the flux executes 
a random walk in time. The emission-line light curves in the 
upper panels exhibit variations correlated with those in the 
driving light curve, but delayed and smeared due to the ap- 
propriate delay maps shown to the left of each emission-line 
light curve. 

A7 Synthetic Spectra 

Fig. A6 shows synthetic spectra calculated for the LOC 
model in the bright and faint states. The velocity distribu- 
tion at e ach radi us is taken to be a Gaussian with dispersion 
At> = ^GM/R for M = 3 x 10 7 M Q . The model spectra bear 
a satisfying resemblance to the observed spectra of Seyfert 1 
galaxies, for example NGC 5548 (Clavel et al. 1991). In par- 
ticular, the C iv/Lya ratio is close to 1. The main differences 
are significantly weaker Mg II and C ill] emission in the pre- 
dicted spectrum. From Fig. A3, we see that the Mgn, Cm], 
and Lya emission regions extend to the 20 light day outer 
radius assumed in this LOC model. These lines would be 
stronger if we increased the outer radius. 

APPENDIX B: THE INVERSE PROBLEM 

Appendix A dealt with the relatively straightforward prob- 
lem of calculating predicted data F\(\,t) given a particu- 



© ???? RAS, MNRAS 000, 1-21 



Quasar Tomography: Unification of Echo Mapping and Photoionisation Models 17 



delay t MEMBLR x /968=1 .01 TEST=0. 00000 

10203040 




KDHl 31-DEC- I'j'jS 15:19 



Figure A5. Synthetic light curves for the LOC model. 



Bright and Faint Spectn: 




Figure A6. Synthetic bright and faint spectra for the LOC 
model. 



(Marsh & Home 1988), and echo mapping (Home, Welsh & 
Peterson 1991, Home 1994). 

Our discussion of the inverse problem is rather formal 
because the maximum entropy fitting techniques are based 
on probability theory and therefore very generally applicable 
rather than being "ad-hoc" methods limited to our specific 
problem. Specific examples of maximum entropy fitting ap- 
plied to the quasar tomography problem are discussed in 
§2. 

We adopt the notation / to refer to a vector (or map) 
of parameters fi, with i labeling the axes of the parameter 
space (or the pixels of the map). In quasar tomography, 
/ includes several positive additive distributions: the cloud 
map f(R, 9, uh, Nk, v), the light curve L(t), the background 
spectrum Fb(X), and the distance D. In practice of course 
we adopt a partition for the domain of each distribution, 
thereby replacing the infinite-dimensional distribution by a 
finite-dimensional map /. 

Constraints on / arise from measurements of a data 
vector D, which impose a probability distribution P(D). In 
quasar tomography, P(D) arises from the incomplete and 
noisy measurements of the time- variable spectrum F\(X,t). 
The forward problem, calculating predicted data R(f) given 
the map /, is well defined (Appendix A). The inverse prob- 
lem, estimating P(f) given P(D), is often under-constrained 
because the number of parameters M = dim(/) exceeds the 
number of data constraints N = dim(Z)). This indetermi- 
nacy can be particularly severe when high-quality datasets 
warrant the use of fine partitions to resolve small structures 
in the maps. 

To cope with the indeterminacy of the inverse prob- 
lem, we regularize the problem by seeking "the simplest" or 
"most probable" solutions that fit the observations. This is a 
good example of "Occam's razor" : when two models succeed 
equally well in accounting for the data, the simpler model is 
more likely to be true. 



Bl To Fit the Data: Maximize Likelihood 

Bayes's theorem tells us how to assess the probability of a 
model with parameters / in the light of a dataset D. The 
posterior probability of / is 



P(f\D) 



P(D\f) Pjf) 
P(D) 



(Bl) 



Here the likelihood, P(D\f), is the probability that D arises 
if we assume that / is true. The prior, P(f), is the prob- 
ability assigned to / before D was obtained. Finally, the 
divisor 



/ 



P(D) = / P(D\f) P(f) df 



(B2) 



lar cloud map f{R, 0, n H , A H , v). We now tackle the more 
treacherous inverse of this problem, fitting an observed 
dataset to recover the underlying cloud map. We accomplish 
this using maximum entropy fitting methods similar to those 
we have developed for other astro-tomography problems, for 
example eclipse mapping (Home 1985), Doppler tomography 



is included to ensure that P(f\D) is a probability density 
normalized to 1 when integrated over all possible maps /. 

The likelihood is easy to evaluate. To be specific about 
the dataset D, consider N independent measurements Di, 
described by Gaussians with standard deviations ct; centred 
on predicted values Ri(f). The likelihood is then 
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Maximum Likelihood Fitting 



Maximum Entropy Fitting 




/ 
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Figure Bl. The ±<r confidence interval and the posterior prob- 
ability distribution are shown for one parameter of a maximum 
likelihood fit. 



- - A - cxp{- X 2 /2} 

^(^i/) = n Pirn) = x Zd ' 



where 



* 2 = E 



Di-Ri(f) 



and the partition function, 



r N 

/ c X p{- X 2 /2}dD = H(2^) 1/2 
•> i=i 



(B3) 



(B4) 



(B5) 



normalizes the likelihood as a probability density on D. 

The maximum likelihood method assumes that P(f) is 
a uniform distribution. The "best fit" model / is the one 
that maximizes the likelihood. This is equivalent (when Oi 
are known) to minimizing x 2 . A model using M parameters 
to fit N data points should leave residuals with N — M 
degrees of freedom. A suc cessful fit sh ould therefore achieve 
xlm/{N - M) w 1 ± y/2/{N -M). The ±a interval for 
any parameter of the fit can be found by the criterion \ 2 < 
Xmin + 1- These concepts are illustrated in Fig. Bl. 

The Bayesian formulation includes the prior P(f) to 
make explicit a fundamental ambiguity involved in the inter- 
pretation of data. Different people may use different priors, 
and thereby arrive at different conclusions from the same 
data. This is illustrated in Fig. B2. 




KDHl 31-DEC- 1998 19:54 



Figure B2. In maximum entropy fitting, we minimize \ 2 — 2aS. 
The cntropic prior P(f\a) oc exp {aS} peaks at the default value 
d and has an rms width (ef/a) 1 / 2 . When the regularization pa- 
rameter a is very small, the prior is so wide that the entropy 
S is ignored and the posterior probability becomes P(f\D,a) oc 
exp {— x 2 /2|. As a increases, the entropy becomes more impor- 
tant. P(f\D,a) then shifts toward the default value and becomes 
narrower. The optimum value of a is determined in Fig. B3. 



B2 For a Simple Map: Maximize Entropy 

The maximum entropy method extends maximum likelihood 
fitting to include models that employ huge numbers of pa- 
rameters. The entropy quantifies the number of parameters 
being used in such fits. Maximum entropy fitting seeks the 
simplest model that fits the data. 

In quasar tomography, our M-dimensional parameter 
vector / (equivalently, M-pixel map) arises from partition- 
ing the domains of several positive additive distributions, 
e.g., f(R, 6, n H , A H , v), L(t), F B (X), and D. If we wish to 
ensure that our maps take on only positive values, and that 
our inferences depend on the underlying distributions rep- 
resented by the maps, and not on specific coordinates or 
partitions used to slice up the domain of /, then the only 
consistent way to do this (Skilling 1989) is to assign relative 
probabilities to the various possible maps / according to an 
entropy of the form 



S(f) =J2 m if* ~ d ^h In (fi/di)} 



(B6) 



Here Wi is a weight proportional to the volume of pixel i, 
and di(f) is the default value of 

The entropy is regarded as measuring the "simplicity" 
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of the map. The simplest map maximizes the entropy, and 
therefore satisfies 

M 

0=J| —In^ + E^d- 1 )^- ^ 

k — 1 

Note that the maximum entropy value, S = 0, occurs when 
/i = d; for every pixel i. The entropy is an information- 
like measure of the distance of / from its default map d(f). 
Maximizing the entropy "steers" / toward d, and that is why 
the di are referred to as default values. We discuss in §B7 
our prescription for chosing the default map d(f) to make 
the entropy express a preference for specific symmetries, for 
example smoothness, or axi-symmetry. 



B3 The Entropic Prior 

The prior probability associated with the entropy (Skilling 
1989) is 



P(f\a)=f[p(Ma) = e -^pl 



(B8) 



This introduces a regularization parameter, a > 0, and the 
corresponding partition function, 



(B9) 



Z s (a) = J P(f\a) df = J exp{aS(/)} df , 

which normalizes the entropic prior to a probability density 
on /. 

To gain some insight into the meaning of this expres- 
sion, expand S(f) in a Taylor series about its peak at / = d, 
and truncate this to obtain the quadratic approximation 



dfidfj 



(BIO) 



The corresponding entropic prior is an M-dimensional Gaus- 
sian, 



P(f\a) « V - >- . 

Z s (a) 

The entropy curvature matrix is given in general by 

a 2 s(/) 

dfidfj 



(Bll) 



ddi 



ddi 



+ V M wt. f f fk l) d2dk fk ddk 9dh 1 



For our present purposes, specialize to the case of a fixed 
default map, ddi/dfj <C 1, for which the entropy curvature 
matrix is diagonal, 



d 2 s(f) 



Wi 



J* 



(B13) 



dfidfj 

The entropic prior is then a product of M Gaussian distri- 
butions 



«<p{-f ££i^(/i-*) a } 



P(J\<*)- y . . 

^s(a) 

with the partition function 
J- V ami 



(B14) 



(B15) 



From this we see that the entropic prior admits values in 
the range /» ~ dj (l ± (QJii,di) -1 ' 2 ). For a 3> (widi) -1 , the 
prior confines attention to a narrow range /< ~ di, while 
for a <JC (widi) -1 , a wider range opens to consideration (see 
Fig. B2). In this way a controls the strength of our prior 
conviction that the default values di are correct. If we wish 
to set a to a different value for each pixel, this can be done 
effectively by changing the pixel weights Wi . 

B4 The Maximum Entropy Trajectory 

Given a dataset D, and a regularization parameter a, the 
posterior probability of our map / is 



P(f\D,a) 



exp{aS- X 2 /2} 



(B16) 



(B17) 



Z Q (a) 

where the appropriate partition function is 
Z Q {a)= J exp{aS- X 2 /2} df . 

As a — > co, the entropy term dominates and the prior be- 
comes so narrow that it over-rules any influence of the data. 
In the opposite limit, a — > 0, the prior becomes very wide, 
and we revert to the maximum likelihood fit minimizing \ 2 - 
Thus a parameterizes a 1-dimensional family of maps f(a), 
the maximum entropy trajectory (Gull & Skilling 1990), con- 
necting the maximum likelihood map (a = 0), and the max- 
imum entropy map (a — » co). 

To explore this family of maximum entropy fits, we use 
the MEMSYS algorithm (Skilling & Bryan 1984) to make 
iterative adjustments to / in order to maximize 



Q = aS - X 2 /2. 



(B18) 



For any given value of a, the iteration proceeds until the gra- 
dients of x 2 and S are found to be parallel to machine preci- 
sion, thus assuring that Q is maximized. We begin where we 
must, with a high value of a, and a simple map that fits the 
data poorly. We then gradually lower a to reduce the influ- 
ence of the entropy in comparison with that of the data. For 
each lower value of a, the map develops additional structure 
in order to improve the fit. Thus x 2 decreases monotonically 
as we pass along the maximum entropy trajectory through 
progressively lower values of a. 



( B12 ) B5 Stopping Criteria 



Which value of a should we choose? How strong is our faith 
that the default values are correct? A justifiable strategy 
rejects high values of a that produce poor fits to the data, 
and low values of a that fit the data too well. Based on this 
motivation, a common practice is to select a value of a that 
achieves x 2 /N = 1. This is the stopping criterion we have 
adopted for all the results shown in this paper. 

The x 2 = N stopping criterion gives a conservative fit. 
It under-fits the data a bit and omits some of the finer-scale 
structure in the map. To see this, consider a maximum like- 
lihood fit with M parameters adjusted to fit N data points. 
The best fit is expected to leave residuals with N — M de- 
grees of freedom, i.e. \ 2 ~ {N - M) ± y / 2(N - M). A fit 
that achieves x 2 ~ N therefore under-fits the data. 

A Bayesian stopping criterion aims to maximize the 
posterior probability of a, given by 
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3ayesian Estimation of a 




og a 
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Figure B3. Bayesian determination of the regularization param- 
eter a for the problem illustrated in Fig. B2. If a is large, the 
entropy pulls the solution to the default value, and the probabil- 
ity is then low because the fit to the data is poor. If a is small, 
the entropy is ignored. The fit is then good, but the probability 
is low because there is a penalty for using a very wide prior. The 
probability peak gives the most likely value of a. 



P{a\D) oc P{D\a) P(a) 

= P(a) JP(D\f) P(f\a) df (B19) 
= P(a) 7 ZQ( 7 a ) , . 

This entails some ambiguity because we are free to choose 
the prior P(a). Since a > 0, a plausible choice is P(a) oc 
a -1 , so that the prior uniform in In a. For an example, see 
Fig. B3. 

We expect P(a\D) to have a maximum at some a = a, 
provided P(a) is not too bizarre. To see this, note that for 
a a the emphasis on entropy maximization makes the 
map very rigid. P{a\D) will then be small because \ 2 i s 
large - very simple maps fail to achieve good fits to the data. 
As we decrease a, we move along the maximum entropy 
trajectory through increasingly flexible maps that develop 
more structure to improve the fit. For a <C a, the map 
is so flexible that x' 2 approaches an asymptotic minimum 
value as the fit becomes as good as it can be. In this limit 
P(a\D) — > because Zs(a) oc a~ M ^ 2 . The entropic prior 
is imposing an appropriate penalty (Occam's razor) on very 
flexible maps that deploy too many degrees of freedom to 
achieve their fit. For a = a we achieve the best compromise 
- fitting the data well enough while deploying a minimum 
number of degrees of freedom in the map. 

How many parameters do we employ in a regularized 



fit? The map may have an enormous number of pixels - in 
fact M may be larger than N - so we can't possibly be using 
all M degrees of freedom. In fact we are not, because the 
structure that develops in the map is spatially correlated. 
Because we are minimizing \ 2 ~ 2aS, it is plausible to iden- 
tify — 2aS as the effective number of parameters used in the 
regularized fit. (Note that — 2aS > because a > and 
S < 0.) By analogy with a maximum likelihood fit, then, 
we should choose a so that \ 2 — N + 2aS. This is only a 
heuristic justification, but the result can be justified more 
rigorously (Gull 1989) when the Bayesian stopping criterion 
is used. 



B6 Quantified Uncertainties 

The uncertainty of the map / after fitting to the data D 
is quantified by its posterior probability distribution. Equa- 
tion (BI6) gives this for each value of a, corresponding to 
different strengths of faith in the default map. 

One can combine these by performing a Bayesian aver- 
age over a to obtain 

PU\B) - f PU\D, a ) P(a)D) da 

In many cases, particularly with large numbers of data 
points and map pixels, the probability peak is narrow 
enough to permit the use of quadratic approximations to 
X 2 and S, yielding an M-dimensional Gaussian probability 
density centred on the most probable map on the maximum 
entropy trajectory at f(a). 

Inferences about other parameters of the fit, such as the 
distance, luminosity, and Ho, may be quantified in the same 
way, 



P(H \D) oc P(D\H ) P(H ) 

= P(H ) JP(D\Ho,f) P(f\a) dfP(a)dc 

= P(go) dfP(a) da 

= P(Ho)fz§^P(*) da. 



(B21) 



B7 Default Maps 

How shall we choose defaults? We want the entropy to 
"steer" the solution toward the "simplest" map that fits the 
data. The default map should therefore be set to whatever 
we consider to be the simplest map. But what do we consider 
be a simple map? 

In general, we compute default values as weighted geo- 
metric means of "nearby" map values: 



Mf) = 0X P \ ^2 Bij ln ^ 



(B22) 



.3 = 1 



Think of the symmetric "blur" matrix Bij as a point-spread 
function, decreasing with the "distance" between pixels i 
and j. It is normalized as a probability distribution 

M 

B v = 1 • ( B23 ) 

3 = 1 



This makes d(f) a "blurred" copy of /. Each pixel in / is 
then pulled toward its neighbors. 
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1/M, then the default map is uniform, weights. We have not attempted to explore this possibility, 
and the maximum entropy fit gives the "most uniform" map but may do so in future work, 
that fits the data. This would be appropriate when the prior 
expectation is that no structure is present in /. In problems 
with incomplete data constraints, a uniform default can al- 
low the data to impress bizarre artifacts on the map, for ex- 
ample the ingress-egress arches seen in eclipse maps (Home 
1985). 

It is usually more appropriate to seek smooth maps by 
using "curvature defaults" that steer each pixel value toward 
the geometric mean of its nearest neighbors, 



di = y/fi+ifi-i . (B24) 

With curvature defaults, the entropy of a 1-dimensional map 
f(x), partitioned at intervals Ax, approaches a weighted in- 
tegral of the squared curvature of the logarithm of the map, 



(Ax) 4 w(x)f(x) ( 9 2 In/ 
dx 2 



(B25) 



With constant w(x), this gives preference to f(x) with Gaus- 
sian peaks (constant curvature) and exponential tails (zero 
curvature). With w(x) oc x 2 , the solution is steered toward 
f(x) that are power-laws in x. 

In dealing with multi-dimensional maps, for example 
f(x) — f(R, 9, hh, -/Vh, v), there are of course many different 
curvatures to consider. We compute the default value for 
pixel x by first calculating geometric means of neighbors 
offset on opposite sides in several chosen directions Ax k , and 
then taking a weighted average of those geometric means: 

, E fc (A fc /Aa- fc ) 2 In ^f(x + Ax k )f(x - Ax k ) 

lnd(z) = 2 .(B26) 

}2 k {A k /Ax k ) 

A large A fc promotes structures oriented along the direction 
Ax k . 

We frequently use the default map to steer our solutions 
toward certain preferred symmetries. For example, to pro- 
mote a spherical geometry, we would compute defaults by 
averaging map values over spherical shells. Point-symmetric 
and axi-symmetric geometries project to maps with front- 
to-back symmetry, i.e. symmetric on reflection through the 
Y axis. To promote point- or axi-symmetric geometries, we 
steer each pixel toward the mean of itself and the pixel with 
opposite sign of cos 9. We also apply constraints on the nor- 
malization of the map, for example covering fraction must 
be less than 100%, by imposing the required constraint on 
the default map. In this way the desired symmetry and nor- 
malization are taken on by the map if this is compatible 
with the data. Otherwise the map comes as close as it can 
to the desired symmetry and normalization while still fitting 
the data. 

Finally, we note that it may be possible by making a 
deft choice of default image to suppress to some extent the 
artefacts that form along the time-delay paraboloids. The 
strategy would be to first blur the map along the time de- 
lay direction, and then mix that in with a negative weight 
to form the final default image. With negative weight the 
default will suppress rather than enhance structure in the 
time-delay direction. This apporach has demonstrated some 
success in suppressing similar artefacts in the eclipse map- 
ping problem (Spruit 1994). A difficulty we see with this 
approach is in deciding how strongly to apply the negative 
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